-
Notifications
You must be signed in to change notification settings - Fork 99
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Legendre, Hermite, Laguerre, Chebyshev polynomials #175
Conversation
Nice! |
I havent looked into those yet but I'd recommend to look at the implementation of python's numpy. I found their Hermite evaluation at github.com/numpy/.... |
@simonbyrne I dont get why the checks fail. Is the build process somehow broken? |
@ludoro You can easily implement Hermite, Laguerre, Legendre, Chebyshev polynomials via recurrence relations. They can all be described by the same recurrence relation and I implemented Hermite and Legendre already. Check out my last commits. |
@simonbyrne I dont plan on implementing any new features for this branch. Let me know what you think of it. I need this branch to start working on implementing spherical harmonics.. |
You need to add Polynomials to the test dependencies in the Project.toml (similar to the Test package: https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/Project.toml) |
Done. |
continuous-integration/drone/pr still fails because it only runs julia 1.0. |
Quote from @MikaelSlevinsky
I liked the formula to reduce the degree and implemented that in the most current commit. Do you know similar formulas for the other polynomials? As to the explicit formula: Usually recurrence formulas should be more efficient (in terms of speed as well as accuracy) but maybe there are parameter ranges where it is better to use the explicit formula? Maybe we could benchmark that at some point.. that could also be done in another issue... |
Thanks for taking this on, @PaulXiCao; it's a big task. Most of the functions in this library are specialized to a set of floating-point types, with some fallbacks for other numbers. In this sense, what makes including orthogonal polynomials challenging is not so much evaluating the recurrences, but that doing so defines the scope of orthogonal polynomials for the entire community. According to Wikipedia, Chebyshev polynomials are also defined for non-natural degree by the formula On the other hand, the O(log n) complexity for the degree doubling formulas is significantly faster than the three-term recurrences when x is a square matrix, but it's not clear when comparing with the analytical formula, which has essentially fixed cost independent of degree. For low degree, the degree doubling seems better, but that changes when the degree gets large: julia> chebyshevt(n, x) = cos(n*acos(x))
chebyshevt (generic function with 1 method)
julia> function chebyshevt(n::Integer, x)
if n == 0
return one(x)
elseif n == 1
return x
elseif iseven(n)
Tn2 = chebyshevt(n÷2, x)
return 2*Tn2*Tn2-one(x)
else
Tn2 = chebyshevt(n÷2, x)
Tn2p1 = chebyshevt(n÷2+1, x)
return 2*Tn2*Tn2p1-x
end
end
chebyshevt (generic function with 2 methods)
julia> x = randn(1000, 1000)/sqrt(2000); x = (x+x')/2;
julia> @time chebyshevt(100, x);
0.321130 seconds (120 allocations: 442.509 MiB, 9.06% gc time)
julia> @time chebyshevt(100.0, x);
1.351499 seconds (84 allocations: 420.002 MiB, 15.79% gc time)
julia> @time chebyshevt(10000, x);
2.151449 seconds (728 allocations: 2.697 GiB, 14.03% gc time)
julia> @time chebyshevt(10000.0, x);
1.510139 seconds (96 allocations: 511.555 MiB, 2.31% gc time)
julia> Another tradeoff is that it appears that the degree-doubling formulas for Chebyshev polynomials have a higher error growth rate than the traditional three-term recurrences. All four types of Chebyshev polynomials, T, U, V, W, would satisfy a degree-doubling-like formula because they're essentially trig functions with a variable transformation. The classical orthogonal polynomials have so-called interior and exterior asymptotic approximations, some of which are mathematically rigorous; see the DLMF https://dlmf.nist.gov/18.15. Implementation of these formulas would naturally come with branching in the code but would make evaluating e.g. Legendre polynomials of extreme degree, say P_{1e24}(0.248), feasible. |
I'm also beginning to wonder whether this is the appropriate package for classical orthogonal polynomials instead of something like JuliaMath/ClassicalOrthogonalPolynomials.jl All classical OPs are of hypergeometric type. So they could be evaluated for argument, parameters, and degree that are |
A package like https://github.com/jverzani/SpecialPolynomials.jl? CC: @jverzani |
That repository should use FastGaussQuadrature.jl for faster algorithms than GLR and FastTransforms.jl for the connection problems. Most of these variations on the same ideas are already unified under one tent: ApproxFun or more generally the JuliaApproximation org. |
Also, it has its own implementation of pFq instead of calling HypergeometricFunctions.jl... The argument for these being in SpecialFunctions.jl is it makes it a “canonical” implementation that can be refined, where putting it in yet another package will mean that they will continue to be reimplemented. |
Yes, agree with all of that. I put a notice in the docs for |
SpecialPolynomials.jl a pretty cool looking package! There's a lot to do to make HypergeometricFunctions.jl more user-friendly, so if you have any use cases, please file issues. The package is more than "just the Maclaurin series," though. As an example, if one defines Chebyshev polynomials by Gauss' hypergeometric function, then julia> using HypergeometricFunctions
julia> chebyshevt(n, x) = _₂F₁(-n, n, 1//2, (1-x)/2)
julia> @time chebyshevt(0.123, 0.456)
0.000003 seconds (5 allocations: 176 bytes)
0.9909056375600692 the actual branch that's called is https://github.com/JuliaMath/HypergeometricFunctions.jl/blob/62ca7c6704cc27c96bbcc76e3573bd780b9ba1e8/src/gauss.jl#L16 which is pretty efficient (further defined here https://github.com/JuliaMath/HypergeometricFunctions.jl/blob/62ca7c6704cc27c96bbcc76e3573bd780b9ba1e8/src/specialfunctions.jl#L81). |
Hi @MikaelSlevinsky that is impressive performance. SpecialPolynomials can evaluate the polynomial in that time, but if you include construction is about 5 times slower :( The only thing I was missing with these implementations was being able to use |
I guess someone needs to decide if this pr is going forward or if it should rather go to a different repo.. |
Let's do a poll: Thumbs up for this repos and Thumbs down for it in a different repos |
I will merge the latest changes from master and resolve the conflicts in the hope that this pr will be merged at some point. Hope dies last ;) |
Codecov Report
@@ Coverage Diff @@
## master #175 +/- ##
==========================================
+ Coverage 88.17% 88.69% +0.52%
==========================================
Files 11 12 +1
Lines 2630 2751 +121
==========================================
+ Hits 2319 2440 +121
Misses 311 311
Flags with carried forward coverage won't be shown. Click here to find out more.
Continue to review full report at Codecov.
|
As I've commented on several of these issues, I tend to feel that orthogonal polynomials belong in a separate package. |
Note I just registered the package: https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl This supports some of these functions, but not all. (It was previously called OrthogonalPolynomialsQuasi.jl) |
It seems like the consensus here is that these belong in a separate package and that @dlfivefifty has registered ClassicalOrthogonalPolynomials.jl. I'm closing this issue for that reason. |
Legendre polynomial is implemented via the Bonnet recursion formula.
Documentation + testing is already inlcuded.
This is in relation to issue #124 .